% Moments

% Variable Profits (i.e., profit excluding any fixed costs):
vaD = piD; % purely domestic
vaM = betastuff * phi.^(sg-1).* sdM.^(g *(1-sg)/(eps-1)); % importer-only
vaX =  betastuff *  phi.^(sg-1) + ...
       betastuff2 * phi.^(theta*(sg-1)) .* fx.^(1-theta); % expoter-only
vaXM=  betastuff *  phi.^(sg-1).* sdXM.^(g *(1-sg)/(eps-1)) + ...
       betastuff2 * phi.^(theta*(sg-1)) .* sdXM.^(theta* g *(1-sg)/(eps-1)) .* fx.^(1-theta); % global
va = (S==1) .* vaD + (S==2) .* vaM + (S==3) .* vaX + (S==4) .* vaXM;

% Material shares:
omega = va/sum(va); % materials and sales are both proportional to variable profits

% Fractions of importers, exporters and importer-exporters
Firms_Domestic  = sum ((S==1))/length(sD); Firms_Importers = sum ((S==2))/length(sD);
Firms_Exporters = sum ((S==3))/length(sD); Firms_ImpExp    = sum ((S==4))/length(sD);
 
% Aggregate import and export Share
Aggsi = sum( (1-sD) .* omega ); AggsX = sum( sX .* omega );

% Standard deviations of log materials, import shares, and export shares:
STD_lnva = std(log(va)); STD_si = (var(1-sD))^(1/2); STD_xi = (var(sX))^(1/2);

% Pairwise correlations between log materials, import shares, and export shares:
corr_sixi = corr(1-sD,sX); corr_lnvasi = corr(log(va),1-sD); corr_lnvaxi = corr(log(va),sX);

% Exports (re-scaled) in domestic labor:  
xX = sg* betastuff2 * phi.^(theta*(sg-1)) .* fx.^(1-theta); % exporter-only
xXM = sg* betastuff2 * phi.^(theta*(sg-1)) .* sdXM.^(theta* g *(1-sg)/(eps-1)) .* fx.^(1-theta); % global
TotalExports = mean( (S==3) .* xX + (S==4) .* xXM ); 

% Imports of materials (use that materials are proportional to variable
% profits):
mM = vaM * g .* (1-sdM) * (sg-1); % importer-only
mXM= vaXM * g .* (1-sdXM) * (sg-1); % global
TotalImports = mean( (S==2) .* mM + (S==4) .* mXM );

% Percentiles of distributions of sales, import and export shares (Table 6)
psales = prctile(log(omega* 200000 / 6856),[10 25 50 75 90 95]); % sales divided by total sales adjusted for the number of firms
psi = prctile((1-sD(sD<1))*100,[10 25 50 75 90 95]); % percentiles of import shares conditional on importing
pxi = prctile(sX(sX>0)*100,    [10 25 50 75 90 95]);

% Conditional distributions (Table 7, first two panels)

pomega = prctile(omega,[25 50 75]); % quartiles of sales
sicond = [mean(1-sD(omega<pomega(1))) mean(1-sD(omega>pomega(1) & omega<=pomega(2))) mean(1-sD(omega>pomega(2) & omega<=pomega(3))) mean(1-sD(omega>pomega(3))) ]*100;
xicond = [mean(sX(omega<pomega(1))) mean(sX(omega>pomega(1) & omega<=pomega(2))) mean(sX(omega>pomega(2) & omega<=pomega(3))) mean(sX(omega>pomega(3)))]*100;

% Average export shares within import share
% quartiles, conditional on global status (Table 7, last panel)
si = 1-sD(sD<1 & sX>0); % import shares conditional on importing-exporting
sx = sX(sD<1 & sX>0); % export shares conditional on importing-exporting
psiXM = prctile(si,[25 50 75]); % quartiles of import shares conditional on importing-exporting
xicondsiXM = [mean(sx(si<=psiXM(1))) mean(sx(si>psiXM(1) & si<=psiXM(2))) mean(sx(si>psiXM(2) & si<=psiXM(3))) mean(sx(si>psiXM(3)))]*100;

% Display targeted moments in model vs data

disp('Targeted moments') 
disp('[1] Aggregate import share')
disp(Aggsi), disp(m(1))
disp('[2] Aggregate export share')
disp(AggsX), disp(m(2))
disp('[3] Fraction of importers')
disp(Firms_Importers+Firms_ImpExp), disp(m(3)+m(5))
disp('[4] Fraction of exporters')
disp(Firms_Exporters+Firms_ImpExp), disp(m(4)+m(5))
disp('[5] Fraction global')
disp(Firms_ImpExp), disp(m(5))
disp('[6] Dispersion in materials')
disp(STD_lnva), disp(m(6))
disp('[7] Dispersion in import shares')
disp(STD_si), disp(m(7))
disp('[8] Dispersion in export shares')
disp(STD_xi), disp(m(8))
disp('[9] Correlation materials-import shares')
disp(corr_lnvasi), disp(m(9))
disp('[10] Correlation materials-export shares')
disp(corr_lnvaxi), disp(m(11))
disp('[11] Correlation import-export shares')
disp(corr_sixi), disp(m(10))